#include "sm2_operation.h"
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

// 2^256 - n = 0x10000000000000000000000008dfc2094de39fad4ac440bf6c62abedd
const sm2_uint256_t SM2_256_MINUS_N = { 0xc62abedd, 0xac440bf6,0xde39fad4, 0x8dfc2094, 0x00000000, 0x00000000, 0x00000000,0x00000001 };

// 256位的1和0
const sm2_uint256_t ZERO = {0};
const sm2_uint256_t ONE  = {1};

/* SM2参数 
*/
const sm2_uint256_t SM2_P = { 0xffffffff, 0xffffffff, 0x00000000, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffffe };
const sm2_uint256_t SM2_A = { 0xfffffffc, 0xffffffff, 0x00000000, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffffe };
const sm2_uint256_t SM2_B = { 0x4D940E93, 0xDDBCBD41, 0x15AB8F92, 0xF39789F5, 0xCF6509A7, 0x4D5A9E4B, 0x9D9F5E34, 0x28E9FA9E };
const sm2_uint256_t SM2_N = { 0x39D54123, 0x53BBF409, 0x21C6052B, 0x7203DF6B, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFE };
const sm2_uint256_t Gx    = { 0x334C74C7, 0x715A4589, 0xF2660BE1, 0x8FE30BBF, 0x6A39C994, 0x5F990446, 0x1F198119, 0x32C4AE2C };
const sm2_uint256_t Gy    = { 0x2139F0A0, 0x02DF32E5, 0xC62A4740, 0xD0A9877C, 0x6B692153, 0x59BDCEE3, 0xF4F6779C, 0xBC3736A2 };

const sm2_uint512_t SM2_P_sqr = { 0x1, 0, 0xfffffffe, 0x1, 0x1, 0xfffffffe, 0, 0x2, 0xfffffffe, 0xfffffffd, 0x3, 0xfffffffe, 0xffffffff, 0xffffffff, 0, 0xfffffffe };   

const sm2_uint256_t SM2_N_MINUS_ONE = { 0x39D54122, 0x53BBF409, 0x21C6052B, 0x7203DF6B, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFE };

// 1P to 13P
const sm2_uint512_t SM2_Ps[13] = {
    {0xffffffff, 0xffffffff, 0, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffffe, 0, 0, 0, 0, 0, 0, 0, 0 },
    {0xfffffffe, 0xffffffff, 0x1, 0xfffffffe, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffffd, 0x1, 0, 0, 0, 0, 0, 0, 0 },
    {0xfffffffd, 0xffffffff, 0x2, 0xfffffffd, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffffc, 0x2, 0, 0, 0, 0, 0, 0, 0 },
    {0xfffffffc, 0xffffffff, 0x3, 0xfffffffc, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffffb, 0x3, 0, 0, 0, 0, 0, 0, 0 },
    {0xfffffffb, 0xffffffff, 0x4, 0xfffffffb, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffffa, 0x4, 0, 0, 0, 0, 0, 0, 0 },
    {0xfffffffa, 0xffffffff, 0x5, 0xfffffffa, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffff9, 0x5, 0, 0, 0, 0, 0, 0, 0 },
    {0xfffffff9, 0xffffffff, 0x6, 0xfffffff9, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffff8, 0x6, 0, 0, 0, 0, 0, 0, 0 },
    {0xfffffff8, 0xffffffff, 0x7, 0xfffffff8, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffff7, 0x7, 0, 0, 0, 0, 0, 0, 0 },
    {0xfffffff7, 0xffffffff, 0x8, 0xfffffff7, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffff6, 0x8, 0, 0, 0, 0, 0, 0, 0 },
    {0xfffffff6, 0xffffffff, 0x9, 0xfffffff6, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffff5, 0x9, 0, 0, 0, 0, 0, 0, 0 },
    {0xfffffff5, 0xffffffff, 0xa, 0xfffffff5, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffff4, 0xa, 0, 0, 0, 0, 0, 0, 0 },
    {0xfffffff4, 0xffffffff, 0xb, 0xfffffff4, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffff3, 0xb, 0, 0, 0, 0, 0, 0, 0 },
    {0xfffffff3, 0xffffffff, 0xc, 0xfffffff3, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffff2, 0xc, 0, 0, 0, 0, 0, 0, 0 }
};


/**
 * @brief   将一个256位数设置为零。
 * 
 * @details 此函数将传入的 `sm2_uint256_t` 类型变量的所有元素都设置为零。
 * 
 * @param[in, out]  num  需要被设置为零的256位数
 * 
 * @return 无
 * 
 * @note    无
 * @warning 使用时需要确保传入的 `num` 已正确分配内存。
 * @pre     无
 * @post    数组中的每个元素都被赋值为零。
 */
void sm2_set_zero(sm2_uint256_t num)
{
    for(int i = 0; i < 8; ++i) num[i] = 0;
}

/**
 * @brief   将一个512位数设置为零。
 * 
 * @details 此函数将传入的 `sm2_uint512_t` 类型变量的所有元素都设置为零。
 * 
 * @param[in, out]  num  需要被设置为零的512位数
 * 
 * @return 无返回值
 * 
 * @note    无
 * @warning 使用时需要确保传入的 `num` 已正确分配内存。
 * @pre     无
 * @post    数组中的每个元素都被赋值为零。
 */
void sm2_set_512_zero(sm2_uint512_t num) 
{
    for(int i = 0; i < 16; ++i) num[i] = 0;
}

/**
 * @brief   将一个256位数设置为1。
 * 
 * @details 此函数将传入的 `sm2_uint256_t` 类型变量的第一个元素设置为1，其余元素设置为零。通常用于初始化为1的数值。
 * 
 * @param[in, out]  num  需要被设置为1的256位数
 * 
 * @return 无返回值
 * 
 * @note    该函数适用于处理 `sm2_uint256_t` 类型的数据结构。
 * @warning 使用时需要确保传入的 `num` 已正确分配内存。
 * @pre     无
 * @post    数组的第一个元素被设置为1，其余元素被设置为零。
 */
void sm2_set_one(sm2_uint256_t num)
{
    num[0] = 1;
    for(int i = 1; i < 8; ++i) num[i] = 0;
}

/**
 * @brief   将一个512位数设置为1。
 * 
 * @details 此函数将传入的 `sm2_uint512_t` 类型变量的第一个元素设置为1，其余元素设置为零。通常用于初始化为1的512位数值。
 * 
 * @param[in, out]  num  需要被设置为1的512位数
 * 
 * @return 无返回值
 * 
 * @note    该函数适用于处理 `sm2_uint512_t` 类型的数据结构。
 * @warning 使用时需要确保传入的 `num` 已正确分配内存。
 * @pre     无
 * @post    数组的第一个元素被设置为1，其余元素被设置为零。
 */
void sm2_set_512_one(sm2_uint512_t num)
{
    num[0] = 1;
    for(int i = 1; i < 16; ++i) num[i] = 0;
}

/**
 * @brief   将一个512位数值转换为256位数值。高位数据被truncate
 * 
 * @details 此函数将传入的 `sm2_uint256_t` 数值的前256位拷贝到 `sm2_uint512_t` 类型的数组 `b` 中，并将 `b` 的高位部分置为零。
 * 
 * @param[in]   a  需要转换的256位数值
 * @param[out]  b  转换后的512位数值
 * 
 * @return 无返回值
 * 
 * @note    该函数假设输入的 `a` 是有效的 `sm2_uint256_t` 数值，且输出的 `b` 必须为有效的 `sm2_uint512_t` 数组。
 * @warning 请确保 `b` 在调用前已经分配足够的内存。
 * @pre     无
 * @post    数组 `b` 的前8个元素被填充为 `a` 的值，其余元素被设置为零。
 */
void sm2_512_to_256(sm2_uint256_t result, const sm2_uint512_t n) {
    sm2_set_512_zero(result);
    for(int i = 0; i < 8; ++i) result[i] = n[i];
}

/**
 * @brief   将一个256位数值转换为512位数值。
 * 
 * @details 此函数将传入的 `sm2_uint256_t` 数值的所有位拷贝到 `sm2_uint512_t` 数组 `b` 的前256位， 
 *          并将 `b` 的其余部分置为零。
 * 
 * @param[in]   a  需要转换的256位数值
 * @param[out]  b  转换后的512位数值
 * 
 * @return 无返回值
 * 
 * @note    该函数假设输入的 `a` 是有效的 `sm2_uint256_t` 数值，且输出的 `b` 必须为有效的 `sm2_uint512_t` 数组。
 * @warning 请确保 `b` 在调用前已经分配足够的内存。
 * @pre     无
 * @post    数组 `b` 的前8个元素被填充为 `a` 的值，其余元素被设置为零。
 */
void sm2_256_to_512(sm2_uint512_t result, const sm2_uint256_t n) {
    sm2_set_512_zero(result);
    for(int i = 0; i < 8; ++i) result[i] = n[i];
}

/**
 * @brief   复制一个512位数值到另一个512位数值。
 * 
 * @details 此函数将源 `sm2_uint512_t` 数值 `a` 的所有元素拷贝到目标 `sm2_uint512_t` 数值 `result` 中。
 * 
 * @param[out]  result  目标512位数值
 * @param[in]   a       源512位数值
 * 
 * @return 无返回值
 * 
 * @note    该函数假设 `a` 和 `result` 都是有效的 `sm2_uint512_t` 数组。
 * @warning 请确保 `result` 在调用前已经分配足够的内存。
 * @pre     无
 * @post    数组 `result` 的元素被填充为 `a` 的对应值。
 */
void sm2_512_copy(sm2_uint512_t result, sm2_uint512_t a) {
    for(int i = 0; i < 16; ++i) result[i] = a[i];
}

/**
 * @brief   复制一个256位数值到另一个256位数值。
 * 
 * @details 此函数将源 `sm2_uint256_t` 数值 `a` 的所有元素拷贝到目标 `sm2_uint256_t` 数值 `result` 中。
 * 
 * @param[out]  result  目标256位数值
 * @param[in]   a       源256位数值
 * 
 * @return 无返回值
 * 
 * @note    该函数假设 `a` 和 `result` 都是有效的 `sm2_uint256_t` 数组。
 * @warning 请确保 `result` 在调用前已经分配足够的内存。
 * @pre     无
 * @post    数组 `result` 的元素被填充为 `a` 的对应值。
 */
void sm2_copy(sm2_uint256_t result, const sm2_uint256_t a) {
    for(int i = 0; i < 8; ++i) result[i] = a[i];
}

/**
 * @brief   比较两个256位数值的大小。
 * 
 * @details 此函数比较两个 `sm2_uint256_t` 类型的数值 `a` 和 `b`。从最高位到最低位逐个元素进行比较， 
 *          若发现某个元素不相等，则根据该元素的大小返回比较结果；若所有元素相等，则返回相等的结果。
 * 
 * @param[in]   a  第一个256位数值
 * @param[in]   b  第二个256位数值
 * 
 * @return 
 *         - `SM2_COMPARE_LT`: 如果 `a` 小于 `b`
 *         - `SM2_COMPARE_GT`: 如果 `a` 大于 `b`
 *         - `SM2_COMPARE_EQ`: 如果 `a` 等于 `b`
 * 
 * @note    函数会逐位比较 `a` 和 `b` 的元素，从高位到低位，若有不相等的元素，则返回对应的比较结果。
 * @warning 确保 `a` 和 `b` 的长度是合法的 `sm2_uint256_t` 数组。
 * @pre     无
 * @post    根据比较结果返回对应的宏值（`SM2_COMPARE_LT`, `SM2_COMPARE_GT`, `SM2_COMPARE_EQ`）。
 */
int sm2_compare(const sm2_uint256_t a, const sm2_uint256_t b)
{
    int i;
    for(i = 7; i >= 0; --i) { // shortcut
        if(a[i] == b[i]) continue;
        else {
            if(a[i] < b[i]) return SM2_COMPARE_LT;
            else return SM2_COMPARE_GT;
        }
    }
    return SM2_COMPARE_EQ;
}

/**
 * @brief   比较两个512位数值的大小。
 * 
 * @details 此函数比较两个 `sm2_uint512_t` 类型的数值 `a` 和 `b`。从最高位到最低位逐个元素进行比较， 
 *          若发现某个元素不相等，则根据该元素的大小返回比较结果；若所有元素相等，则返回相等的结果。
 * 
 * @param[in]   a  第一个512位数值
 * @param[in]   b  第二个512位数值
 * 
 * @return 
 *         - `SM2_COMPARE_LT`: 如果 `a` 小于 `b`
 *         - `SM2_COMPARE_GT`: 如果 `a` 大于 `b`
 *         - `SM2_COMPARE_EQ`: 如果 `a` 等于 `b`
 * 
 * @note    函数会逐位比较 `a` 和 `b` 的元素，从高位到低位，若有不相等的元素，则返回对应的比较结果。
 * @warning 确保 `a` 和 `b` 的长度是合法的 `sm2_uint512_t` 数组。
 * @pre     无
 * @post    根据比较结果返回对应的宏值（`SM2_COMPARE_LT`, `SM2_COMPARE_GT`, `SM2_COMPARE_EQ`）。
 */
int sm2_compare_512(const sm2_uint512_t a, const sm2_uint512_t b)
{
    int i;
    for(i = 15; i >= 0; --i) { // shortcut
        if(a[i] == b[i]) continue;
        else {
            if(a[i] < b[i]) return SM2_COMPARE_LT;
            else return SM2_COMPARE_GT;
        }
    }
    return SM2_COMPARE_EQ;
}


/**
 * @brief   执行两个256位数值的加法运算。
 * 
 * @details 此函数计算两个 `sm2_uint256_t` 类型的数值 `a` 和 `b` 的和，并将结果存储在 `result` 中。加法过程
 *          考虑了进位，并逐位进行加法计算。函数返回加法结果中的最终进位值。
 * 
 * @param[out]  result  存储加法结果的 `sm2_uint256_t` 数值
 * @param[in]   a       第一个参与加法运算的256位数值
 * @param[in]   b       第二个参与加法运算的256位数值
 * 
 * @return 
 *         - 返回加法运算后的最终进位值（0或1）
 * 
 * @note    该函数假设输入的 `a`、`b` 和 `result` 都是有效的 `sm2_uint256_t` 数组。
 * @warning 请确保 `result` 在调用前已分配足够的内存空间。
 * @pre     无
 * @post    `result` 数组存储了 `a` 和 `b` 加法运算后的结果。
 */
uint32_t sm2_add(sm2_uint256_t result, const sm2_uint256_t a, const sm2_uint256_t b)
{
    uint32_t t = 0, c = 0;
    sm2_set_zero(result);
    t = a[0] + b[0];
    c = t < a[0];
    result[0] = t;

    for(int i = 1; i < 8; ++i) {
        t = a[i] + c;
        c = t < a[i];
        result[i] = t + b[i];
        c += result[i] < t; //进位
    }

    return c;
}

/**
 * @brief   执行两个512位数值的加法运算。
 * 
 * @details 此函数计算两个 `sm2_uint512_t` 类型的数值 `a` 和 `b` 的和，并将结果存储在 `result` 中。加法过程
 *          考虑了进位，并逐位进行加法计算。函数返回加法结果中的最终进位值。
 * 
 * @param[out]  result  存储加法结果的 `sm2_uint512_t` 数值
 * @param[in]   a       第一个参与加法运算的512位数值
 * @param[in]   b       第二个参与加法运算的512位数值
 * 
 * @return 
 *         - 返回加法运算后的最终进位值（0或1）
 * 
 * @note    该函数假设输入的 `a`、`b` 和 `result` 都是有效的 `sm2_uint512_t` 数组。
 * @warning 请确保 `result` 在调用前已分配足够的内存空间。
 * @pre     无
 * @post    `result` 数组存储了 `a` 和 `b` 加法运算后的结果。
 */
uint32_t sm2_512_add(sm2_uint512_t result, const sm2_uint512_t a, const sm2_uint512_t b)
{
    uint32_t t = 0, c = 0;
    sm2_set_512_zero(result);
    t = a[0] + b[0];
    c = t < a[0];
    result[0] = t;

    for(int i = 1; i < 16; ++i) {
        t = a[i] + c;
        c = t < a[i];
        result[i] = t + b[i];
        c += result[i] < t; // 进位
    }

    return c;
}

/**
 * @brief   执行两个256位数值的减法运算。
 * 
 * @details 此函数计算两个 `sm2_uint256_t` 类型的数值 `a` 和 `b` 的差，并将结果存储在 `result` 中。减法过程
 *          考虑了借位，并逐位进行减法计算。函数返回减法结果中的最终借位值。
 * 
 * @param[out]  result  存储减法结果的 `sm2_uint256_t` 数值
 * @param[in]   a       第一个参与减法运算的256位数值
 * @param[in]   b       第二个参与减法运算的256位数值
 * 
 * @return 
 *         - 返回减法运算后的最终借位值（0或1）
 * 
 * @note    该函数假设输入的 `a`、`b` 和 `result` 都是有效的 `sm2_uint256_t` 数组。
 * @warning 请确保 `result` 在调用前已分配足够的内存空间。
 * @pre     无
 * @post    `result` 数组存储了 `a` 和 `b` 减法运算后的结果。
 */
uint32_t sm2_sub(sm2_uint256_t result, const sm2_uint256_t a, const sm2_uint256_t b) 
{
    sm2_set_zero(result);
    uint32_t t = 0, c = 0;
    t = a[0] - b[0];
    c = t > a[0];
    result[0] = t;

    for(int i = 1; i < 8; ++i) {
        t = a[i] - c;
        c = t > a[i];
        result[i] = t - b[i];
        c += result[i] > t; // 借位
    }

    return c;
}

/**
 * @brief   执行两个512位数值的减法运算。
 * 
 * @details 此函数计算两个 `sm2_uint512_t` 类型的数值 `a` 和 `b` 的差，并将结果存储在 `result` 中。减法过程
 *          考虑了借位，并逐位进行减法计算。函数返回减法结果中的最终借位值。
 * 
 * @param[out]  result  存储减法结果的 `sm2_uint512_t` 数值
 * @param[in]   a       第一个参与减法运算的512位数值
 * @param[in]   b       第二个参与减法运算的512位数值
 * 
 * @return 
 *         - 返回减法运算后的最终借位值（0或1）
 * 
 * @note    该函数假设输入的 `a`、`b` 和 `result` 都是有效的 `sm2_uint512_t` 数组。
 * @warning 请确保 `result` 在调用前已分配足够的内存空间。
 * @pre     无
 * @post    `result` 数组存储了 `a` 和 `b` 减法运算后的结果。
 */
uint32_t sm2_512_sub(sm2_uint512_t result, const sm2_uint512_t a, const sm2_uint512_t b) {
    uint32_t t = 0, c = 0;
    t = a[0] - b[0];
    c = t > a[0];
    result[0] = t;

    for(int i = 1; i < 16; ++i) {
        t = a[i] - c;
        c = t > a[i];
        result[i] = t - b[i];
        c += (uint32_t)(result[i] > t);
    }

    return c;
}

/**
 * @brief   执行两个256位数值的乘法运算，结果存储在512位数值 `r` 中。
 * 
 * @details 此函数通过逐位计算，使用16位单位对两个 `sm2_uint256_t` 类型的数值 `a` 和 `b` 进行乘法运算。
 *          计算结果存储在 `r` 中。乘法过程分为以下几步：
 *          1. 将 `a` 和 `b` 拆分成低16位和高16位两个部分。
 *          2. 对每个16位部分进行逐步乘法运算。
 *          3. 结果累加至数组 `s` 中，并考虑进位。
 *          4. 最终将结果拼接并赋值给输出 `r`。
 * 
 * @param[out]  r       存储乘法结果的 `sm2_uint512_t` 数值
 * @param[in]   a       第一个参与乘法运算的256位数值
 * @param[in]   b       第二个参与乘法运算的256位数值
 * 
 * @return 无返回值
 * 
 * @note    该函数使用了16位的单位进行运算，因此输入的 `a` 和 `b` 会被拆分成更小的部分进行处理。
 * @warning 请确保 `r` 在调用前已分配足够的内存空间，并且 `a` 和 `b` 都是有效的 `sm2_uint256_t` 数组。
 * @pre     无
 * @post    `r` 数组存储了 `a` 和 `b` 乘法运算的结果。
 */
void sm2_mul(sm2_uint512_t r, const sm2_uint256_t a, const sm2_uint256_t b)
{
    uint32_t a_[16] = {0}; // 16bit unit 数组用于存储 `a` 的低16位和高16位部分
    uint32_t b_[16] = {0}; // 16bit unit 数组用于存储 `b` 的低16位和高16位部分
    uint32_t c;           // 用于存储乘法运算中的进位
    uint32_t s[32] = {0}; // 存储部分结果的数组，最终的乘法结果将累加到这里
    sm2_set_512_zero(r);  // 初始化结果数组 r，确保其初始值为 0

    // 将 a 和 b 拆分为低16位和高16位
    for(int i = 0; i < 8; ++i) {
        a_[i*2] = a[i] & 0xffff;     // 存储 a 的低 16 位部分
        b_[i*2] = b[i] & 0xffff;     // 存储 b 的低 16 位部分
        a_[i*2+1] = a[i] >> 16;      // 存储 a 的高 16 位部分
        b_[i*2+1] = b[i] >> 16;      // 存储 b 的高 16 位部分
    }

    // 进行逐位的乘法运算，并累加到结果数组 s 中
    for(int i = 0; i < 16; ++i) {
        c = 0;  // 每次开始新一轮乘法时，初始化进位为 0
        for(int j = 0; j < 16; ++j) {
            // 计算当前位的乘法结果并累加到 s[i+j] 中，考虑进位
            c = s[i+j] + a_[i] * b_[j] + c;
            s[i+j] = c & 0xffff;
            c >>= 16;
        }
        s[i+16] = c;  // 处理进位并存储
    }

    // 将最终的结果拼接到 r 中
    for(int i = 0; i < 16; ++i) {
        r[i] |= (s[i*2] & 0xffff);       // 获取 s 中的低 16 位并存储到 r 中
        r[i] |= (s[i*2+1] << 16);        // 获取 s 中的高 16 位并存储到 r 中
    }
}

/**
 * @brief   执行两个256位数值的加法并进行模运算。
 * 
 * @details 此函数首先计算两个 `sm2_uint256_t` 类型的数值 `a` 和 `b` 的和，并将结果存储在 `result` 中。 
 *          然后，结果通过模运算（使用 `sm2_fast_mod_reduction` 函数）进行规范化，确保最终的结果符合规定的范围。
 * 
 * @param[out]  result  存储加法并模运算后的结果的 `sm2_uint256_t` 数值
 * @param[in]   a       第一个参与加法运算的256位数值
 * @param[in]   b       第二个参与加法运算的256位数值
 * 
 * @return 无返回值
 * 
 * @note    该函数假设 `a`、`b` 和 `result` 都是有效的 `sm2_uint256_t` 数组。
 * @warning 请确保 `result` 在调用前已分配足够的内存空间。
 * @pre     无
 * @post    `result` 数组存储了 `a` 和 `b` 加法运算并经过模运算后的结果。
 */
void sm2_mod_add(sm2_uint256_t result, const sm2_uint256_t a, const sm2_uint256_t b)
{
    sm2_uint512_t tmp = {0};
    int c = sm2_add(tmp, a, b);
    tmp[8] = c; // 保留进位
    sm2_fast_mod_reduction(result, tmp); // 快速模P约减
}


/**
 * @brief   执行两个256位数值的减法运算，并进行模运算。
 * 
 * @details 此函数首先比较两个 `sm2_uint256_t` 类型的数值 `a` 和 `b`，若 `a` 小于 `b`，则先用模数 `SM2_P` 
 *          对 `b` 进行减法运算，并将结果与 `a` 相加；否则直接进行 `a` 和 `b` 的减法运算。最后将结果存储在 `result` 中。
 *          函数返回比较结果，若 `a < b`，返回 `SM2_COMPARE_LT`，否则执行减法并返回无进位。
 * 
 * @param[out]  result  存储减法并模运算后的结果的 `sm2_uint256_t` 数值
 * @param[in]   a       第一个参与减法运算的256位数值
 * @param[in]   b       第二个参与减法运算的256位数值
 * 
 * @return 
 *         - `SM2_COMPARE_LT`: 如果 `a` 小于 `b`，则返回此值
 *         - 无返回值：直接执行减法并将结果存储在 `result` 中
 * 
 * @note    该函数使用了预定义的常量 `SM2_P`，确保它是合法的并已定义。
 * @warning 请确保 `result` 在调用前已分配足够的内存空间，并且 `a` 和 `b` 都是有效的 `sm2_uint256_t` 数组。
 * @pre     无
 * @post    `result` 数组存储了 `a` 和 `b` 减法运算并经过模运算后的结果。
 */
void sm2_mod_sub(sm2_uint256_t result, const sm2_uint256_t a, const sm2_uint256_t b)
{
    sm2_uint256_t tmp = {0};
    if(sm2_compare(a, b) == SM2_COMPARE_LT){
        sm2_sub(tmp, SM2_P, b);
        sm2_add(result, a, tmp);
        return;
    }
    sm2_sub(result, a, b);
}

/**
 * @brief   计算给定256位数值的模负值。
 * 
 * @details 此函数计算给定的 `sm2_uint256_t` 类型数值 `a` 的模负值，具体通过将模数 `SM2_P` 减去 `a` 来实现，
 *          并将结果存储在 `result` 中。即，`result = SM2_P - a`。
 * 
 * @param[out]  result  存储模负值结果的 `sm2_uint256_t` 数值
 * @param[in]   a       参与计算模负值的256位数值
 * 
 * @return 无返回值
 * 
 * @note    该函数使用了预定义的常量 `SM2_P`，确保它是合法的并已定义。
 * @warning 请确保 `result` 在调用前已分配足够的内存空间，并且 `a` 是有效的 `sm2_uint256_t` 数组。
 * @pre     无
 * @post    `result` 数组存储了 `a` 的模负值。
 */
void sm2_mod_neg(sm2_uint256_t result, const sm2_uint256_t a)
{
    sm2_sub(result, SM2_P, a);      // result = -a mod P = P - a
}

/**
 * @brief   计算两个256位数值 `a` 和 `b` 的模乘积，并将结果存储在 `r` 中。
 * 
 * @details 此函数首先计算 `a` 和 `b` 的乘积，并将结果存储在 `tmp` 中。然后通过快速模约简将结果
 *          映射到一个256位的数值，并将最终结果存储在 `r` 中。该过程保证了返回的结果是小于一个
 *          固定的模数（如椭圆曲线的模数）。
 * 
 * @param[out]  r       存储模乘法结果的 `sm2_uint256_t` 数值
 * @param[in]   a       第一个参与乘法的256位数值
 * @param[in]   b       第二个参与乘法的256位数值
 * 
 * @return 无返回值
 * 
 * @note    调用该函数时，结果 `r` 必须为一个预先分配的有效数组，并且 `a` 和 `b` 应为有效的
 *          `sm2_uint256_t` 数组。
 * @pre     `a` 和 `b` 都是有效的 `sm2_uint256_t` 数组
 * @post    `r` 数组存储了 `a` 和 `b` 模乘的结果。
 */
void sm2_mod_mul(sm2_uint256_t r, const sm2_uint256_t a, const sm2_uint256_t b)
{
    sm2_uint512_t tmp = {0};
    sm2_mul(tmp, a, b);
    sm2_fast_mod_reduction(r, tmp);
}

/**
 * @brief   计算一个256位数值 `a` 的模平方，并将结果存储在 `result` 中。
 * 
 * @details 此函数通过调用 `sm2_mod_mul` 函数计算 `a` 和 `a` 的模乘积（即模平方）。结果会被存储
 *          在 `result` 中。通过这种方式，可以有效地计算平方值并保持结果在有效范围内。
 * 
 * @param[out]  result  存储模平方结果的 `sm2_uint256_t` 数值
 * @param[in]   a       需要进行模平方的256位数值
 * 
 * @return 无返回值
 * 
 * @note    `a` 和 `result` 都应为有效的 `sm2_uint256_t` 数组。
 * @pre     `a` 是一个有效的 `sm2_uint256_t` 数组
 * @post    `result` 数组存储了 `a` 的模平方结果。
 */
void sm2_mod_sqr(sm2_uint256_t result, const sm2_uint256_t a) 
{
    sm2_mod_mul(result, a, a);
}


/**
 * @brief   计算32bit常数 `k` 与256位数值 `a` 的模乘积，并将结果存储在 `result` 中。
 * 
 * @details 此函数首先将常数 `k` 转换为一个256位数值 `_k`，然后计算 `_k` 与 `a` 的模乘积。
 *          结果通过调用 `sm2_mod_mul` 存储在 `result` 中。该过程确保结果是模约简后的正确值。
 * 
 * @param[out]  result  存储模乘法结果的 `sm2_uint256_t` 数值
 * @param[in]   k       常数 `k`，将与 `a` 进行模乘法运算
 * @param[in]   a       参与运算的256位数值
 * 
 * @return 无返回值
 * 
 * @note    `result` 和 `a` 应为有效的 `sm2_uint256_t` 数组，`k` 是32位常数。
 * @pre     `a` 和 `result` 都是有效的 `sm2_uint256_t` 数组
 * @post    `result` 数组存储了常数 `k` 和 `a` 的模乘积结果。
 */
void sm2_mod_k_mul(sm2_uint256_t result, const uint32_t k, const sm2_uint256_t a)
{
    sm2_uint256_t _k = {0};
    _k[0] = k;
    sm2_mod_mul(result, _k, a);
}

/**
 * @brief   快速模约简运算，使用 SM2 算法的特定素数参数。
 * 
 * @details 该函数通过将一个512位的大整数 `a` 按照 SM2 算法的模约简方法进行约简，最终将结果存储在256位的 `result` 中。
 *          运算过程涉及到对一个大整数 `a` 进行多次相加、相减操作，并通过与 SM2 模数的比较来调整最终结果。
 *          参考白国强等人的《SM2椭圆曲线公钥密码算法》中的模约简过程。
 *          Z. Zhao and G. Bai, "Ultra High-Speed SM2 ASIC Implementation," 2014 IEEE 13th International Conference on Trust, Security and Privacy in Computing and Communications, Beijing, China, 2014, pp. 182-188, doi: 10.1109/TrustCom.2014.27.
 *          keywords: {Computer architecture;Elliptic curve cryptography;Elliptic curves;Pipelines;Jacobian matrices;Algorithm design and analysis;Optimization;Elliptic Curve Cryptography;SM2;Point Multiplication;High Speed;ASIC},
 * 
 * @param[out] result  存储模约简结果的 `sm2_uint256_t` 数值
 * @param[in]  a       参与模约简的512位大整数
 * 
 * @return 无返回值
 * 
 * @note    `result` 和 `a` 应为有效的 `sm2_uint256_t` 或 `sm2_uint512_t` 数组，`a` 为512位大整数。
 * @pre     `a` 是有效的 `sm2_uint512_t` 数组。
 * @post    `result` 数组存储了 `a` 模约简后的256位结果。
 */
void sm2_fast_mod_reduction(sm2_uint256_t result, const sm2_uint512_t a)
{
    // SM2素数参数 Psm2p-256 = 2^256-2^224-2^96+2^64-1
    // c是一个大整数(512bit) c = (c[15], c[14],..., c[1], c[0]), (0 <= c < p^2)
    int i = 0;

    sm2_uint512_t _result = {0};  // 临时存储约简结果

    // 定义14个256位数值，按照SM2的模约简过程进行分解
    sm2_uint256_t s[14] = {
        {a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7]},
        {a[8], a[9], 0, a[11], a[12], a[13], a[14], a[15]},
        {a[13], a[14], 0, a[13], a[14], a[15], 0, a[14]},
        {a[14], a[15], 0, 0, 0, 0, 0, a[13]},
        {a[15], 0, 0, 0, 0, 0, 0, a[12]},
        {a[12], a[13], 0, a[14], a[15], a[10], a[11], a[11]},
        {a[10], a[11], 0, a[12], a[13], a[14], a[15], a[10]},
        {a[9], a[10], 0, a[8], a[9], 0, 0, a[9]},
        {a[11], a[12], 0, a[15], 0, 0, 0, a[8]},
        {0, 0, 0, 0, 0, 0, 0, a[15]},
        {0, 0, a[14], 0, 0, 0, 0, 0},
        {0, 0, a[13], 0, 0, 0, 0, 0},
        {0, 0, a[9], 0, 0, 0, 0, 0},
        {0, 0, a[8], 0, 0, 0, 0, 0},
    };
    
    uint32_t c = 0; // 临时变量用于存储进位
    sm2_uint256_t sum = {0};  // 存储相加结果
    sm2_uint512_t _sum = {0}; // 存储最终结果
    sm2_uint256_t temp = {0}; // 临时存储每次运算的中间结果

    // 依次加法相加，每次计算加法并更新 sum
    for(i = 0; i < 10; ++i) {
        c += sm2_add(temp, sum, s[i]);  // 加法并更新进位
        sm2_copy(sum, temp);  // 将加法结果保存到 sum
    }

    c += sm2_add(temp, sum, s[2]);
    sm2_copy(sum, temp);
    c += sm2_add(temp, sum, s[3]);
    sm2_copy(sum, temp);
    c += sm2_add(temp, sum, s[4]);
    sm2_copy(sum, temp);
    c += sm2_add(temp, sum, s[9]);
    sm2_copy(sum, temp);
    c -= sm2_sub(temp, sum, s[10]);
    sm2_copy(sum, temp);
    c -= sm2_sub(temp, sum, s[11]);
    sm2_copy(sum, temp);
    c -= sm2_sub(temp, sum, s[12]);
    sm2_copy(sum, temp);
    c -= sm2_sub(temp, sum, s[13]);
    sm2_copy(sum, temp);

    sm2_256_to_512(_sum, sum);          // 将 sum 从256位转换为512位
    _sum[8] = c;                        // 补上 sum 运算的最终进位

    // 比较 _sum 和 SM2 的模数，进行模约简
    if(sm2_compare_512(_sum, SM2_Ps[0]) == SM2_COMPARE_LT) sm2_512_copy(_result, _sum);  // 如果 _sum 小于第一个模数，直接复制
    else {
        // 否则，通过比较找到合适的模数，并进行减法操作
        for(i = 1; i < 13; ++i) {
            if(sm2_compare_512(_sum, SM2_Ps[i]) == SM2_COMPARE_LT) {
                break;
            }
        }
        sm2_512_sub(_result, _sum, SM2_Ps[i-1]);  // 减去合适的模数，模数的范围：0, 1P, 2P, ..., 13P
    }

    // 将最终结果的低256位存储到result中
    for(int i = 0; i < 8; ++i) result[i] = _result[i];
}

/**
 * @brief   计算 SM2 的模逆元。
 * 
 * @details 该函数通过多次平方、乘法运算以及模约简运算，使用 SM2 算法中的特定方法计算输入 `a` 的模逆元，并将结果存储在 `result` 中。
 *          运算过程包括对 `a` 的多次幂运算，以及通过连续的平方和乘法步骤得到最终的逆元。
 * 
 * @param[out] result  存储计算结果的 `sm2_uint256_t` 数值
 * @param[in]  a       需要计算模逆的 `sm2_uint256_t` 数值
 * 
 * @return 无返回值
 * 
 * @note    `result` 和 `a` 应为有效的 `sm2_uint256_t` 数组，`a` 为256位数值。
 * @pre     `a` 是有效的 `sm2_uint256_t` 数组。
 * @post    `result` 数组存储了 `a` 的模逆元。
 */
void sm2_fast_mod_inverse(sm2_uint256_t result, const sm2_uint256_t a)
{
    sm2_uint256_t temp1={0}, temp2={0}, temp3={0}, temp4={0};
    sm2_uint256_t tmp={0};
    sm2_uint256_t big_tmp={0};
    int i;

    sm2_set_zero(result);
    sm2_mod_sqr(tmp, a);
    sm2_mod_mul(temp1, a, tmp);         // temp1 = a^2 * a

    sm2_mod_sqr(tmp, temp1);
    sm2_mod_mul(temp2, tmp, a);         // temp2 = temp1^2 * a

    sm2_copy(tmp, temp2);
    for(i=0; i<3; ++i){
        sm2_mod_sqr(temp3, tmp);
        sm2_copy(tmp, temp3);
    }
    sm2_mod_mul(temp3, tmp, temp2);     // temp3 = temp2 ^ (2^3) * temp2

    sm2_mod_sqr(temp2, temp3);          // temp2 = temp3 ^ 2

    sm2_mod_mul(temp4, temp2, a);       // temp4 = temp2 * a

    sm2_copy(tmp, temp2);
    for (i = 0; i < 5; i++) {
        sm2_mod_sqr(temp2, tmp);
        sm2_copy(tmp, temp2);
    }          
    sm2_mod_mul(temp2, tmp, temp3);     // temp2 =  temp2 ^ (2^5) * temp3

    sm2_copy(tmp, temp2);
    for (i = 0; i < 12; i++) {
        sm2_mod_sqr(temp3, tmp);
        sm2_copy(tmp, temp3);
    }          
    sm2_mod_mul(temp3, tmp, temp2);     // temp3 = temp2 ^ (2^12) * temp2

    sm2_copy(tmp, temp3);
    for (i = 0; i < 7; i++) {
        sm2_mod_sqr(temp3, tmp);
        sm2_copy(tmp, temp3);
        
    }          
    sm2_mod_mul(temp2, tmp, temp4);     // temp2 = temp3 ^ (2^7) * temp4 

    sm2_mod_sqr(tmp, temp2);
    sm2_mod_sqr(temp3, tmp);            // temp3 = temp2 ^ (2^2)

    sm2_copy(temp4, temp3);
    for (i = 0; i < 29; i++) {
        sm2_mod_sqr(tmp, temp4);
        sm2_copy(temp4, tmp);
    }                                   // temp4 = temp3 ^ (2^29)

    sm2_copy(tmp, temp2);
    sm2_mod_mul(temp2, tmp, temp4);     // temp2 = temp2 * temp4

    sm2_copy(tmp, temp4);
    sm2_mod_sqr(temp4, tmp);
    sm2_copy(tmp, temp4);
    sm2_mod_sqr(temp4, tmp);            // temp4 = temp4 ^ (2^2)

    sm2_mod_mul(tmp, temp1, temp4);
    sm2_mod_mul(temp1, temp3, tmp);     // temp1 = temp1 * temp4 * temp3
    
    sm2_copy(tmp, temp4);
    for (i = 0; i < 32; i++) {
        sm2_mod_sqr(big_tmp, tmp);
        sm2_copy(tmp, big_tmp);
    }          
    sm2_copy(big_tmp, tmp);
    sm2_mod_mul(tmp, big_tmp, temp1);
    for (i = 0; i < 64; i++) {
        sm2_copy(big_tmp, tmp);
        sm2_mod_sqr(tmp, big_tmp);
    }
    sm2_copy(big_tmp, tmp);
    sm2_mod_mul(tmp, big_tmp, temp1);
    for (i = 0; i < 94; i++) {
        sm2_copy(big_tmp, tmp);
        sm2_mod_sqr(tmp, big_tmp);
    }
    sm2_copy(temp3, tmp);               // temp3 = ( (temp4 ^ (2^32) * temp1 )^(2^64) * temp1  )^(2^94)
    
    sm2_mod_mul(tmp, temp2, temp3);
    sm2_copy(big_tmp, tmp);
    sm2_mod_sqr(tmp, big_tmp);
    sm2_copy(big_tmp, tmp);
    sm2_mod_sqr(tmp, big_tmp);
    sm2_mod_mul(result, tmp, a);        // result = temp1 = (temp2 * temp3)^(2^2) * a
}

/**
 * @brief   计算 SM2 算法中的模 N 加法。
 * 
 * @details 该函数实现了模 N 加法运算，即计算 `(a + b) mod N`，其中 `N` 为 SM2 算法中的常量。
 *          如果 `a + b` 大于 N，则减去 N。否则，如果 `a + b` 小于 2^256，则加上 `2^256 - N`。
 * 
 * @param[out] result  存储计算结果的 `sm2_uint256_t` 数值
 * @param[in]  a       加数 `a`
 * @param[in]  b       加数 `b`
 * 
 * @return 无返回值
 * 
 * @note    `result` 和 `a`、`b` 应为有效的 `sm2_uint256_t` 数组，且其长度为 256 位。
 * @pre     `a` 和 `b` 是有效的 `sm2_uint256_t` 数组。
 * @post    `result` 数组存储了 `(a + b) mod N` 的结果。
 */
void sm2_modn_add(sm2_uint256_t result, const sm2_uint256_t a, const sm2_uint256_t b)
{
    int c = sm2_add(result, a, b);
    sm2_uint256_t tmp;
    if(c){
        // a + b - n = (a + b - 2^256) + (2^256 - n)
        sm2_add(tmp, result, SM2_256_MINUS_N);
        sm2_copy(result, tmp);
    }else if(sm2_compare(result, SM2_N) == SM2_COMPARE_GT){
        sm2_sub(tmp, result, SM2_N);
        sm2_copy(result, tmp);
    }
}

/**
 * @brief   计算 SM2 算法中的模 N 减法。
 * 
 * @details 该函数实现了模 N 减法运算，即计算 `(a - b) mod N`，其中 `N` 为 SM2 算法中的常量。
 *          如果 `a - b` 小于 0，则加上 `2^256 - N`。否则，如果 `a - b` 大于 N，则减去 N。
 * 
 * @param[out] result  存储计算结果的 `sm2_uint256_t` 数值
 * @param[in]  a       被减数 `a`
 * @param[in]  b       减数 `b`
 * 
 * @return 无返回值
 * 
 * @note    `result` 和 `a`、`b` 应为有效的 `sm2_uint256_t` 数组，且其长度为 256 位。
 * @pre     `a` 和 `b` 是有效的 `sm2_uint256_t` 数组。
 * @post    `result` 数组存储了 `(a - b) mod N` 的结果。
 */
void sm2_modn_sub(sm2_uint256_t result, const sm2_uint256_t a, const sm2_uint256_t b)
{
    int c = sm2_sub(result, a, b);
    sm2_uint256_t tmp;
    if(c){
        // a - b + n = (a - b + 2^256) - (2^256 - n)
        sm2_sub(tmp, result, SM2_256_MINUS_N);
        sm2_copy(result, tmp);
    }else if(sm2_compare(result, SM2_N) == SM2_COMPARE_GT){
        sm2_sub(tmp, result, SM2_N);
        sm2_copy(result, tmp);
    }
}


/**
 * @brief   计算 SM2 算法中的模 N 取负运算。
 * 
 * @details 该函数实现了模 N 取负运算，即计算 `(-a) mod N`，等价于 `SM2_N - a`。
 * 
 * @param[out] result  存储计算结果的 `sm2_uint256_t` 数值
 * @param[in]  a       输入数值 `a`
 * 
 * @return 无返回值
 * 
 * @note    `result` 和 `a` 应为有效的 `sm2_uint256_t` 数组，且其长度为 256 位。
 * @pre     `a` 是有效的 `sm2_uint256_t` 数组。
 * @post    `result` 数组存储了 `(-a) mod N` 的结果。
 */
void sm2_modn_neg(sm2_uint256_t result, const sm2_uint256_t a)
{
    sm2_modn_sub(result, SM2_N, a);
}


// n' = -n^(-1) mod 2^256
//    = 0x6f39132f82e4c7bc2b0068d3b08941d4df1e8d34fc8319a5327f9e8872350975
// sage: -(IntegerModRing(2^256)(n))^-1

const sm2_uint256_t SM2_NEG_N_INV_MOD_Z256 = {
    0x72350975,0x327f9e88,0xfc8319a5,0xdf1e8d34,0xb08941d4,0x2b0068d3,0x82e4c7bc,0x6f39132f
};

// mont(1) (mod n) = 2^256 - n
const uint32_t *SM2_MODN_MONT_ONE = SM2_256_MINUS_N;


/**
 * @brief   执行 SM2 算法中的蒙特卡罗乘法（Montgomery multiplication）模 N 运算。
 * 
 * @details 该函数实现了在模 N 下的蒙特卡罗乘法，即计算 `(a * b) mod N`，采用蒙特卡罗方法优化运算。
 * 
 * @param[out] result  存储计算结果的 `sm2_uint256_t` 数值
 * @param[in]  a       输入数值 `a`
 * @param[in]  b       输入数值 `b`
 * 
 * @return 无返回值
 * 
 * @note    `result` 和 `a`, `b` 应为有效的 `sm2_uint256_t` 数组，且其长度为 256 位。
 * @pre     `a` 和 `b` 是有效的 `sm2_uint256_t` 数组。
 * @post    `result` 数组存储了 `(a * b) mod N` 的结果。
 */
void sm2_modn_mont_mul(sm2_uint256_t result, const sm2_uint256_t a, const sm2_uint256_t b)
{
    sm2_uint512_t z, t;
    uint32_t c;
    sm2_uint256_t tmp;
    sm2_uint512_t tmp1;

    sm2_mul(z, a, b); // z = a * b
    sm2_mul(t, z, SM2_NEG_N_INV_MOD_Z256); // t = low(z) * (-n^(-1) mod 2^256)
    sm2_copy(tmp, t);
    sm2_mul(t, tmp, SM2_N); // t = low(t) * n
    sm2_512_copy(tmp1, z);
    c = sm2_512_add(z, tmp1, t); // z = z + t
    sm2_copy(result, z+8); // result = hign(result)
    if(c){
        sm2_add(tmp, result, SM2_MODN_MONT_ONE);
        sm2_copy(result, tmp);
    }else if (sm2_compare(result, SM2_N) == SM2_COMPARE_GT){
        sm2_sub(tmp, result, SM2_N);
        sm2_copy(result, tmp);
    }
}

/**
 * @brief   执行 SM2 算法中的模 N 乘法。
 * 
 * @details 该函数实现了在模 N 下的乘法运算 `(a * b) mod N`，通过先将输入数值转换为蒙特卡罗域，再进行蒙特卡罗乘法操作，最后将结果转换回常规域。
 * 
 * @param[out] result  存储计算结果的 `sm2_uint256_t` 数值
 * @param[in]  a       输入数值 `a`
 * @param[in]  b       输入数值 `b`
 * 
 * @return 无返回值
 * 
 * @note    `result`, `a`, `b` 应为有效的 `sm2_uint256_t` 数组，且其长度为 256 位。
 * @pre     `a` 和 `b` 是有效的 `sm2_uint256_t` 数组。
 * @post    `result` 数组存储了 `(a * b) mod N` 的结果。
 */
void sm2_modn_mul(sm2_uint256_t result, const sm2_uint256_t a, const sm2_uint256_t b)
{
	sm2_uint256_t _a;
	sm2_uint256_t _b;

	sm2_modn_to_mont(a, _a);
	sm2_modn_to_mont(b, _b);
	sm2_modn_mont_mul(result, _a, _b);
	sm2_modn_from_mont(result, result);
}


/**
 * @brief   执行 SM2 算法中的模 N 蒙特卡罗平方运算。
 * 
 * @details 该函数实现了在模 N 下的平方运算 `(a * a) mod N`，通过调用模 N 蒙特卡罗乘法函数来实现平方操作。
 * 
 * @param[out] result  存储计算结果的 `sm2_uint256_t` 数值
 * @param[in]  a       输入数值 `a`
 * 
 * @return 无返回值
 * 
 * @note    `result`, `a` 应为有效的 `sm2_uint256_t` 数组，且其长度为 256 位。
 * @pre     `a` 是有效的 `sm2_uint256_t` 数组。
 * @post    `result` 数组存储了 `(a * a) mod N` 的结果。
 */
void sm2_modn_mont_sqr(sm2_uint256_t result, const sm2_uint256_t a)
{
    sm2_uint256_t tmp = {0};
    sm2_copy(tmp, a);
	sm2_modn_mont_mul(result, tmp, a);
}

/**
 * @brief   执行 SM2 算法中的模 N 平方运算。
 * 
 * @details 该函数先将输入数值转换为蒙特卡罗域，再执行平方操作，最后将结果转换回常规域。
 * 
 * @param[out] result  存储计算结果的 `sm2_uint256_t` 数值
 * @param[in]  a       输入数值 `a`
 * 
 * @return 无返回值
 * 
 * @note    `result`, `a` 应为有效的 `sm2_uint256_t` 数组，且其长度为 256 位。
 * @pre     `a` 是有效的 `sm2_uint256_t` 数组。
 * @post    `result` 数组存储了 `(a * a) mod N` 的结果。
 */
void sm2_modn_sqr(sm2_uint256_t result, const sm2_uint256_t a)
{
	sm2_uint256_t _a;

	sm2_modn_to_mont(a, _a);
	sm2_modn_mont_sqr(result, _a);
	sm2_modn_from_mont(result, result);
}

/**
 * @brief   执行 SM2 算法中的模 N 蒙特卡罗指数运算。
 * 
 * @details 该函数实现了模 N 的蒙特卡罗指数运算 `a^b mod N`。它首先将常数值 `1` 转换为蒙特卡罗域形式，并依次计算 `a` 的各个幂次，通过平方和乘法运算来完成指数运算。
 * 
 * @param[out] result  存储计算结果的 `sm2_uint256_t` 数值
 * @param[in]  a       输入基数 `a`
 * @param[in]  b       输入指数 `b`
 * 
 * @return 无返回值
 * 
 * @note    `result`, `a`, `b` 应为有效的 `sm2_uint256_t` 数组，且其长度为 256 位。
 * @pre     `a`, `b` 是有效的 `sm2_uint256_t` 数组，`b` 表示指数。
 * @post    `result` 数组存储了 `a^b mod N` 的结果。
 */
void sm2_modn_mont_exp(sm2_uint256_t result, const sm2_uint256_t a, const sm2_uint256_t b)
{
	sm2_uint256_t t = {0};
    sm2_uint256_t _t = {0};
	uint32_t w[2];
	int i, j;

	// t = mont(1)
	sm2_copy(t, SM2_MODN_MONT_ONE);

	for (i = 3; i >= 0; i--) {
		// w = b[i];
        w[0] = b[i<<1];
        w[1] = b[(i<<1)+1];
		for (j = 0; j < 64; j++) {
            sm2_copy(_t, t);
			sm2_modn_mont_sqr(t, _t);
			if (w[1] & 0x80000000) {
                sm2_copy(_t, t);
				sm2_modn_mont_mul(t, _t, a);
			}
            // w <<= 1;
			w[1] <<= 1;
            w[1] |= w[0] >> 31;
            w[0] <<= 1;
		}
	}

	sm2_copy(result, t);
}

/**
 * @brief   执行 SM2 算法中的模 N 指数运算。
 * 
 * @details 该函数实现了模 N 指数运算 `a^b mod N`。首先将基数 `a` 转换为蒙特卡罗域形式，然后调用 `sm2_modn_mont_exp` 执行模 N 的蒙特卡罗指数运算，最后将结果从蒙特卡罗域形式转换回常规形式。
 * 
 * @param[out] result  存储计算结果的 `sm2_uint256_t` 数值
 * @param[in]  a       输入基数 `a`
 * @param[in]  b       输入指数 `b`
 * 
 * @return 无返回值
 * 
 * @note    `result`, `a`, `b` 应为有效的 `sm2_uint256_t` 数组，且其长度为 256 位。
 * @pre     `a`, `b` 是有效的 `sm2_uint256_t` 数组，`b` 表示指数。
 * @post    `result` 数组存储了 `a^b mod N` 的结果。
 */
void sm2_modn_exp(sm2_uint256_t result, const sm2_uint256_t a, const sm2_uint256_t b)
{
	sm2_uint256_t _a = {0};

	sm2_modn_to_mont(a, _a);
	sm2_modn_mont_exp(result, _a, b);
    sm2_copy(_a, result);
	sm2_modn_from_mont(result, _a);
}

// n - 2 = 0xfffffffeffffffffffffffffffffffff7203df6b21c6052b53bbf40939d54121
const sm2_uint256_t SM2_N_MINUS_TWO = {
	0x39d54121, 0x53bbf409, 0x21c6052b, 0x7203df6b, 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffffe
};
// TODO: use the special form of SM2_Z256_N_MINUS_TWO[2, 3]		

/**
 * @brief   执行 SM2 算法中的模 N 蒙特卡罗反元素运算。
 * 
 * @details 该函数实现了计算 SM2 算法中的模 N 蒙特卡罗反元素运算，即计算给定数值 `a` 的模 N 逆元。通过反复平方和乘法运算，使用蒙特卡罗方法计算反元素并将结果存储在 `result` 中。
 * 
 * @param[out] result  存储计算结果的 `sm2_uint256_t` 数值
 * @param[in]  a       输入数值 `a`，用于计算其模 N 逆元
 * 
 * @return 无返回值
 * 
 * @note    `result`, `a` 应为有效的 `sm2_uint256_t` 数组，且其长度为 256 位。
 * @pre     `a` 是有效的 `sm2_uint256_t` 数组，表示待求逆元的数值。
 * @post    `result` 数组存储了 `a` 在模 N 下的逆元。
 */
void sm2_modn_mont_inverse(sm2_uint256_t result, const sm2_uint256_t a)
{
	sm2_uint256_t t;
	uint32_t w[2];
	int i;
	int k = 0;

	sm2_copy(t, a);

	for (i = 0; i < 30; i++) {
		sm2_modn_mont_sqr(t, t);
		sm2_modn_mont_mul(t, t, a);
	}
	sm2_modn_mont_sqr(t, t);
	for (i = 0; i < 96; i++) {
		sm2_modn_mont_sqr(t, t);
		sm2_modn_mont_mul(t, t, a);
	}
	w[0] = SM2_N_MINUS_TWO[2];
    w[1] = SM2_N_MINUS_TWO[3];
	for (i = 0; i < 64; i++) {
		sm2_modn_mont_sqr(t, t);
		if (w[1] & 0x80000000) {
			sm2_modn_mont_mul(t, t, a);
		}
		// w <<= 1;
        w[1] <<= 1;
        w[1] |= w[0] >> 31;
        w[0] <<= 1;
	}
	w[0] = SM2_N_MINUS_TWO[0];
    w[1] = SM2_N_MINUS_TWO[1];
	for (i = 0; i < 64; i++) {
		sm2_modn_mont_sqr(t, t);
		if (w[1] & 0x80000000) {
			sm2_modn_mont_mul(t, t, a);
		}
		// w <<= 1;
        w[1] <<= 1;
        w[1] |= w[0] >> 31;
        w[0] <<= 1;
	}

	sm2_copy(result, t);
}

/**
 * @brief   执行 SM2 算法中的模 N 逆元素运算。
 * 
 * @details 该函数实现了计算 SM2 算法中的模 N 逆元素运算。通过将输入数值 `a` 转换为蒙特卡罗表示，调用 `sm2_modn_mont_inverse` 计算其模 N 逆元，并将结果转换回常规表示存储在 `result` 中。
 * 
 * @param[out] result  存储计算结果的 `sm2_uint256_t` 数值
 * @param[in]  a       输入数值 `a`，用于计算其模 N 逆元
 * 
 * @return 无返回值
 * 
 * @note    `result`, `a` 应为有效的 `sm2_uint256_t` 数组，且其长度为 256 位。
 * @pre     `a` 是有效的 `sm2_uint256_t` 数组，表示待求逆元的数值。
 * @post    `result` 数组存储了 `a` 在模 N 下的逆元。
 */
void sm2_modn_inverse(sm2_uint256_t result, const sm2_uint256_t a)
{
	sm2_uint256_t _a;

	sm2_modn_to_mont(a, _a);
	sm2_modn_mont_inverse(result, _a);
	sm2_modn_from_mont(result, result);
}


// mont(mont(a), 1) = aR * 1 * R^-1 (mod n) = a (mod p)
/**
 * @brief   将数值从蒙特卡罗表示转换回常规表示。
 * 
 * @details 该函数通过将输入数值 `a` 与常量 `ONE` 相乘，执行模 N 下的蒙特卡罗转换。计算结果存储在 `result` 中。此操作用于从蒙特卡罗表示恢复原始数值。
 * 
 * @param[out] result  存储转换结果的 `sm2_uint256_t` 数值
 * @param[in]  a       输入数值 `a`，需要从蒙特卡罗表示转换的数值
 * 
 * @return 无返回值
 * 
 * @note    `result`, `a` 应为有效的 `sm2_uint256_t` 数组，且其长度为 256 位。
 * @pre     `a` 是有效的 `sm2_uint256_t` 数组，表示蒙特卡罗表示的数值。
 * @post    `result` 数组存储了转换回常规表示后的数值。
 */
void sm2_modn_from_mont(sm2_uint256_t result, const sm2_uint256_t a)
{
	sm2_modn_mont_mul(result, a, ONE);
}

// 2^512 (mod n) = 0x1eb5e412a22b3d3b620fc84c3affe0d43464504ade6fa2fa901192af7c114f20
const sm2_uint256_t SM2_512_MOD_N = {
	0x7c114f20, 0x901192af, 0xde6fa2fa, 0x3464504a, 0x3affe0d4, 0x620fc84c, 0xa22b3d3b, 0x1eb5e412
};

/**
 * @brief   将数值从常规表示转换为蒙特卡罗表示。
 * 
 * @details 该函数通过将输入数值 `a` 与常量 `SM2_512_MOD_N` 相乘，执行模 N 下的蒙特卡罗转换。计算结果存储在 `r` 中。此操作用于将常规表示的数值转换为蒙特卡罗表示。
 * 
 * @param[in]  a       输入数值 `a`，需要转换为蒙特卡罗表示的数值
 * @param[out] r       存储转换结果的 `sm2_uint256_t` 数值
 * 
 * @return 无返回值
 * 
 * @note    `a`, `r` 应为有效的 `sm2_uint256_t` 数组，且其长度为 256 位。
 * @pre     `a` 是有效的 `sm2_uint256_t` 数组，表示常规表示的数值。
 * @post    `r` 数组存储了转换后的蒙特卡罗表示数值。
 */
void sm2_modn_to_mont(const sm2_uint256_t a, sm2_uint256_t r)
{
    // mont(a) = a * 2^256 (mod n) = mont_mul(a, 2^512 mod n)
	sm2_modn_mont_mul(r, a, SM2_512_MOD_N);
}


// 一个工具函数，按照小端法打印sm2_uint256_t的值
void print_256(const sm2_uint256_t n)
{
    int i = 0;
    for(i = 0; i < 7; ++i) {
        printf("%#.8x, ", n[i]);
    }
    printf("%#.8x", n[i]);
    printf("\n");
}

// 一个工具函数，按照小端法打印sm2_uint512_t的值
void print_512(const sm2_uint512_t n)
{
    int i = 0;
    for(i = 0; i < 15; ++i) {
        printf("%#.8x, ", n[i]);
    }
    printf("%#.8x", n[i]);
    printf("\n");
}

int sm2_print(FILE *fp, int ind, int fmt, const char *label, const sm2_uint256_t a)
{
    fprintf(fp, "%*s%s: ", ind, "", label);
    for (int i = 7; i >= 0; i--)
    {
        fprintf(fp, "%08x", a[i]);
    }
    fprintf(fp, "\n");
    return 1;
}

/*
int main(void)
{
    //sm2_uint512_t num1_0 = {0x700dcbf6, 0x87a68461, 0xdb3d92fe, 0x194a4832, 0x5f3cfcb1, 0xf7a57dae, 0xc65dc547, 0xbf3402f8, 0xc78387ed, 0xd524bcf0, 0x0fbb123f, 0x48b12505, 0x47f8bc7b, 0x74f5cdcc, 0x7b49bb6f, 0xcead1358};
    
    //sm2_uint512_t num1_0 = {0xc78387ed, 0xd524bcf0, 0x0fbb123f, 0x48b12505, 0x47f8bc7b,  0x7b49bb6f, 0xcead1358};
    sm2_uint512_t r;
    sm2_uint512_t one;
    sm2_set_512_one(one);
    sm2_uint512_t result;
    sm2_512_add(result, SM2_Ps[2], one);
    sm2_fast_mod_reduction(r, result);

}
*/

/*
int main(void)
{
    sm2_uint256_t Ps = {0};
    Ps[0] = 13;
    sm2_uint512_t result;

    printf("P^2: \n");
    sm2_mul(result, Ps, SM2_P);
    printf("{");
    for(int j = 0; j < 15; ++j) {
        printf("%#x, ", result[j]);
    }
    printf("%#x ", result[15]);
    printf("};\n");
}
*/
    /* For testing*/
    /*
    sm2_uint256_t num1_0 = {0x700dcbf6, 0x87a68461, 0xdb3d92fe, 0x194a4832, 0x5f3cfcb1, 0xf7a57dae, 0xc65dc547, 0xbf3402f8};
    sm2_uint256_t num2_0 = {0xc78387ed, 0xd524bcf0, 0x0fbb123f, 0x48b12505, 0x47f8bc7b, 0x74f5cdcc, 0x7b49bb6f, 0xcead1358};
    
    sm2_uint256_t num1_1 = {0xf05fc1b3, 0xb9b33a69, 0x9fc2f9db, 0x525230a0, 0x92afb669, 0x2d7f59d2, 0xcce04c9d, 0x1e910287};
    sm2_uint256_t num2_1 = {0x7956656b, 0x481ef69a, 0x783cb3de, 0xd8822204, 0xde7a254e, 0xcbbe9bd2, 0x7b1e05e3, 0x533984fb};
    
    sm2_uint256_t num1_2 = {0x82f78887, 0xb80d6886, 0x23807711, 0x352cedb1, 0x93a2ee46, 0x216d1323, 0xd422e6f9, 0x16ae83e0};
    sm2_uint256_t num2_2 = {0xbeec1906, 0x962887c5, 0x3769f04c, 0x631d6aed, 0xcc985c95, 0xfa3f00cf, 0xde4553ed, 0xaa8157a2};
    
    sm2_uint256_t num1_3 = {0xc297a500, 0x727ca28e, 0x41234711, 0x8c04d313, 0x202d3ab6, 0x91c78bee, 0x6f3c0f91, 0xd6ff0e82};
    sm2_uint256_t num2_3 = {0xdd85875b, 0x7ccb142a, 0x336102b5, 0x4c143bda, 0xc071e3af, 0x1a167cb0, 0x9e91922d, 0x806af7d3};
    
    sm2_uint256_t num1_4 = {0xa577a5cb, 0x10b46b13, 0x74e55594, 0xd93937c8, 0x94051338, 0x08993f06, 0x31311ee5, 0xb2c18a8f};
    sm2_uint256_t num2_4 = {0x335685f4, 0x8979a0ec, 0xd6aba79f, 0x05ee9000, 0xf0c78e99, 0x3b9013b6, 0x78238b3f, 0x7edb07e2};
    
    sm2_uint256_t num1_5 = {0x5add5109, 0x643434c7, 0x1bde90aa, 0x0da9ea73, 0xdbbe8ecb, 0xab2ac40f, 0x77f2f3f3, 0x45a7d5fd};
    sm2_uint256_t num2_5 = {0x36d2e633, 0x729a946a, 0x549ae7b5, 0xc5dc4be9, 0x40f77df8, 0xa017269b, 0x86b1fde8, 0x4ba7e96a};
    
    sm2_uint256_t num1_6 = {0xcfbf1036, 0xeceb3388, 0x74f44fdf, 0xe04bbcb2, 0x3ae43949, 0x107cb020, 0x1b564421, 0xf396ac6e};
    sm2_uint256_t num2_6 = {0x279cba13, 0xf8bf38b3, 0x33a2ca50, 0xcca0dd25, 0xf17d42ae, 0x936fce4a, 0xc7419ce7, 0x271459cb};
    
    sm2_uint256_t num1_7 = {0xf5421b66, 0x582ce6e1, 0xe40582bb, 0x8f2d8438, 0x88f06074, 0x452c3a4f, 0x3bd0cf16, 0x7b481d6f};
    sm2_uint256_t num2_7 = {0xfc8e91bd, 0x41739e04, 0x596801a2, 0x73cba247, 0xfbc7e2e9, 0x101a7233, 0xb797abe2, 0x92d729b2};
    
    sm2_uint256_t num1_8 = {0xa45b566b, 0xc959aa07, 0x933e8d10, 0xe7c2db39, 0x0c5ac441, 0x85ec8ea1, 0x40ac2dee, 0x16502826};
    sm2_uint256_t num2_8 = {0x25c7d6a8, 0x97f0ac4d, 0x69505664, 0x16a41144, 0x5913e2a8, 0x729a8e60, 0x0566250f, 0x79767622};
    
    sm2_uint256_t num1_9 = {0x98b80721, 0xc7fa87a3, 0x9ad33d00, 0x1ce4c5d7, 0xf6a7b1d3, 0x585f852b, 0xfb55f38a, 0x9ef7cb7d};
    sm2_uint256_t num2_9 = {0x61a8b3f3, 0x4ec36d17, 0x81ba812d, 0x0ed36bfe, 0xfeba4a33, 0x8e734dea, 0x5c2cb91f, 0x00bef935};

    sm2_uint512_t result;
    sm2_mul(result, num1_0, num2_0);
    for(int i = 0; i < 16; ++i) {
        printf("%#x, ", result[i]);
    }
    printf("\n");
    */

